rm(list = ls())
options(stringsAsFactors = F)
library(tidyverse)
library(clusterProfiler)
library(msigdbr)  #install.packages("msigdbr")
library(GSVA) 
library(GSEABase)
library(pheatmap)
library(limma)
library(BiocParallel)
## msigdbr包提取下载 一般尝试KEGG和GO做GSVA分析
##KEGG
KEGG_df_all <-  msigdbr(species = "Homo sapiens", # Homo sapiens or Mus musculus
                        category = "C2",
                        subcategory = "CP:KEGG") 
KEGG_df <- dplyr::select(KEGG_df_all,gs_name,gs_exact_source,gene_symbol)
kegg_list <- split(KEGG_df$gene_symbol, KEGG_df$gs_name) ##按照gs_name给gene_symbol分组

##GO
GO_df_all <- msigdbr(species = "Homo sapiens",
                     category = "C5")  
GO_df <- dplyr::select(GO_df_all, gs_name, gene_symbol, gs_exact_source, gs_subcat)
GO_df <- GO_df[GO_df$gs_subcat!="HPO",]
go_list <- split(GO_df$gene_symbol, GO_df$gs_name) ##按照gs_name给gene_symbol分组

##提取自己的样本
FPKM=read.table("convert_exp.txt",sep="\t",header=T,check.names=F, row.names = 1)
####  GSVA  ####
#GSVA算法需要处理logCPM, logRPKM,logTPM数据或counts数据的矩阵####
#dat <- as.matrix(counts)
#dat <- as.matrix(log2(edgeR::cpm(counts))+1)
#dat <- as.matrix(log2(tpm+1))
dat <- as.matrix(log2(FPKM+1))
##这一步是提取自己的基因集
genelist <- read.table("工作簿3.csv",sep=",",header=T,check.names=F)
geneset <- lapply(genelist, function(col) col[!is.na(col) & col != ""])

#geneset <- go_list##这一步是用go的

##linux系统
gsva_mat <- gsva(expr=dat, 
                 gset.idx.list=geneset, 
                 kcdf="Gaussian" ,#"Gaussian" for logCPM,logRPKM,logTPM, "Poisson" for counts
                 verbose=T, 
                 parallel.sz = parallel::detectCores())#调用所有核
##win系统

#param <- SnowParam(workers = parallel::detectCores(), type = "SOCK")
##gsva_mat <- gsva(expr = dat, gset.idx.list = geneset, kcdf = "Gaussian", BPPARAM = param)
##CPU会直接占满卡死，不能用
counts=read.table("GSE126044_counts.txt",sep="\t",header=T,check.names=F, row.names = 1)
####  GSVA  ####
#GSVA算法需要处理logCPM, logRPKM,logTPM数据或counts数据的矩阵####
dat <- as.matrix(counts)
#dat <- as.matrix(log2(edgeR::cpm(counts))+1)
#dat <- as.matrix(log2(tpm+1))
#dat <- as.matrix(log2(FPKM+1))
##这一步是提取自己的基因集
genelist <- read.table("工作簿3.csv",sep=",",header=T,check.names=F)
geneset <- lapply(genelist, function(col) col[!is.na(col) & col != ""])

#geneset <- go_list##这一步是用go的

##linux系统
gsva_mat2 <- gsva(expr=dat, 
                 gset.idx.list=geneset, 
                 kcdf="Gaussian" ,#"Gaussian" for logCPM,logRPKM,logTPM, "Poisson" for counts
                 verbose=T, 
                 parallel.sz = parallel::detectCores())#调用所有核
##win系统
merge <- cbind(gsva_mat2,gsva_mat)
write.csv(gsva_mat,"gsva_go_matrix.csv")
